Statistical investigations on nitrogen-vacancy center creation 
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Experiments show that shallow nitrogen implantations (<10keV) result in a negatively charged 
nitrogen- vacancy center (NV - ) yield of 0.01 — 0.1%. The most succesful technique for introducing 
NV - centers in the carbon matrix is ion implantation followed by annealing at HOOK. We investi- 
gated the influence of channeling effects during shallow implantation and statistical diffusion using 
molecular dynamics (MD) and Monte Carlo (MC) approaches. Energy barriers for the diffusion 
process were calculated using the density functional theory (DFT). Our simulations show a signif- 
icant difference in the NV yield compared to the experiment. Statistically, 25% of the implanted 
nitrogens form a NV center after annealing. 
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I. INTRODUCTION 

A basic requirement for quantum information pro- 
cessing (QIP), magnetometry and even for biomarkers 
is the accurate placement of impurities in the target 
material. A promising candidate for these applications 
is the negatively charged nitrogen- vacancy (NV - ) 
center embedded in diamond. This defect consists of a 
carbon (C) atom which is replaced by a nitrogen (N) 
together with an adjacent vacancy (V) next to it and 
an additional electron captured from the diamond lattice. 

Experimentaly, ion implantation techniques are used to 
introduce N and V defects in the carbon matrix. Upon 
annealing the diamond at HOOK 1 for several hours, 
the fixed N captures mobile V creating NV centers. 
Experiments show that shallow nitrogen implantations 
(<10keV) result in a NV" yield of 0.01 - 0.1%2. 

Here we present simulation results of NV center 
creation. Initially, we studied the statistical influence of 
channeling effects^ on shallow implantations. In further 
investigations, we focused on the statistical diffusion of 
vacancies in diamond and its impact on the NV center 
yield. Our simulations show a significant difference in 
the NV yield compared to experiments. Simulations 
suggest that 25% of the implanted nitrogens build a NV 
center after annealing. 



For atom distances r < 0.3 A the Tersoff potential is 
smoothly switched to the ZBL potential. Functions 
for turning-on and turning-off potentials are given by 
f on = 1 - e - 12 ^+ - 05 ) 2 and f off = e - 12 (-+ - 05 ) 2 where r 
is the distance of the atoms. A Nose-Hoover thermostat 
was used to simulate the temperature^—. 

In order to bring a system into equilibrium at room 
temperature, the thermal cycle shown in figure [1] was 
employed: 
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Figure 1: Schematical initialization plot: Temperature versus 
time including six simulation steps 



II. METHOD 

Using the ITAP Molecular Dynamics (IMD) 
Programme a diamond system of 5.5 x 5.5 x 52nm 3 
was produced. The C-C and N-C interaction were rep- 
resented by using a Tersoff potential^— with interaction 
parameters from refs^^. For decreasing C-C or N-C 
distances the potential energy of the Tersoff potential 
reaches a constant limiting value allowing the atoms to 
move arbitarily close to one other. Therefore, a Ziegler- 
Biersack-Littmark (ZBL) potential was adde d 11 ! 12 . 



It consists of five NPT and one final NVE ensemble 
simulation steps. The first four are used to heat the 
structure up to 1000K, relax it at 1000K, cool it down 
to room temperature (RT) and to relax the diamond at 
RT using periodic boundaries in all directions. In the 
last NPT step the periodic boundary conditions in the 
z direction (implantation direction) are removed. The 
final NVE simulation equilibrates the system by holding 
the ground layer of the diamond fixed in the z direction. 
Figure [T] shows a plot of the diamond temperature 
over the simulated time. The six simulation steps are 
indicated in the plot. 
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All NPT ensemble simulations have been done with a 
time step of 0.1 fs and the NVE equilibration with time 
steps of O.Olfs. The implantation process itself was 
simulated by adding the impact nitrogen atom on top 
of the diamond surface ((100) direction), at a distance 
d = 2.51A above the surface; this distance is beyond 
the range of the Tersoff potential. The particle is then 
set upon a normal incident trajectory with the kinetic 
energy Ek = 4keV. 

We analyse the structure after implantation by com- 
paring the damaged structure with the initial structure 
before implantation (reference structure). Around each 
atom position of the reference structure a shell with a 
radius of r = 1.54A is placed. Counting the number 
of atoms of the damaged structure in current shell 
allows to analyze the integrity of the diamond lattice 
(0 atoms= Vacancy; l=Lattice atom; 2=Interstitial 
atom and lattice atom; 2 > Complex structure). We 
also examine the impact atom related information: 
Depth of the impact atom, displacement in ^-direction 
(straggling), kinetic energy of the impact atom. 

In addition Crystal-TRIM (CTRIM)^ calculations 
are compared to MD results. The CTRIM method uses 
a ZBL potential for the atom interactions and includes 
an electronic stopping power which is missing in the MD 
calculations. 

To investigate the probability for NV center cre- 
ation during annealing, we also conducted kinetic 
Monte-Carlo simulations based on the simple hopping 
frequency of defects: 

T = uj a e~^. (1) 

Here co a = 10 13 s _1 denotes the attempt frequency with 
which defects attempt to overcome a barrier of energy 
E a . The values for those activation energies were taken 
from the refs^ii"— . 

Additionally, DFT simulations using the ABINIT— ~— 
package have been conducted to investigate a possible 
uphill potential a vacancy has to climb when approach- 
ing a nitrogen atom. Such a potential could consider- 
ably decrease the possibility of NV center formation. For 
the calculations the projector augmented wave (PAW) 
method and a 2 x 2 x 2 Monkhorst-Pack grid were used 
for k-point sampling. The exchange correlation poten- 
tial (V xc ) was approximated by local density approxima- 
tion (LDA). The examined structures, including second, 
third, fourth and sixth next neighbor configurations of 
an on-site nitrogen and a vacancy, showed the same total 
energy within errors. The only exception is a specific con- 
figuration with a vacancy on the most distant third next 
neighbor positions (figure [2]). With a nitrogen on lattice 
site A it is energetically more favorable for a vacancy to 
be at the third next-neigbor site labelled 1 instead of 2. 
Here, due to symmetry reasons the lattice is able to relax 



which decreases the total energy by about 0.7eV. There- 
fore no hint for such an uphill potential could be found. 




Figure 2: Diamond structure with a nitrogen impurity (A) 
and a vacancy labeled by 1 or 2. The total energy of the anti- 
symmetric configuration is 0.7eV lower than of the symmetric 
configuration. 



We considered defect structures of different implanta- 
tion energies. Since IMD does not take the electronic 
stopping into account only the 4keV defect structure is 
simulated by this method. For higher implantation en- 
ergies defect structures simulated by the Stopping and 
Range of Ions in Matter (SRIM)±±i2i (5-50keV) are used. 
The crystal size chosen in the simulations corresponds to 
an implantation fluence of 1 • 10 12 cm -2 dose which would 
produce a concentration of about 1 • 10 17 cm -3 nitrogen 
atoms. A simulated time of 2 hours of annealing at HOOK 
decreases the dynamics of vacancies within the crystal in 
a way that only a negligible number of them is still un- 
dergoing diffusion at the end of the simulation. 



III. RESULTS 

A. Implantation 

In recent years, surface impurities, in particular NV 
centers in diamond, have been extensively studied exper- 
imentally. Magnetic field imaging 25 , fluorescence reso- 
nance energy transfer (FRET)2£ and nuclear magnetic 
resonance (NMR) 27 using NV - centers in diamond rely 
on the weak magnetic dipole interaction which decreases 
with the distance by 1/r 3 . Accordingly, applications in- 
volving the fine structure splitting in NV - require shal- 
low implantation, so that the NV - experiences as large 
a field as possible. 

Several studies showed a significant influence of chan- 
neling effects on the depth distribution of implanted ni- 
trogens into diamon d 24 ! 28 . Usual methods as SRIM un- 
derestimate the defect depth by a factor of two 28 . The 
reason can be seen in the fact that SRIM assumes an 
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amorphous structure for the diamond. We study implan- 
tation processes using IMD and CTRIM which take the 
crystal structure into account. 

Figure [3] shows the depth distribution of 4keV nitrogen 
implantations into diamond at 0° with respect to the axis 
of incidence for (100). The very time demanding MD cal- 
culations were repeated 120 times to provide statistics. 
All nitrogen atoms have a different initial lateral start- 
ing point in order to get statistical average values for 
the depth distribution. In comparison adequate CTRIM 
simulations were done. 



the crystal, i.e. in CTRIM one implants into the pre 
damaged lattice. Additionaly the empirical potential 
used in both methods give rise to different results. The 
bond ordered Tersoff potential used in MD describes 
the interaction between the projectile and the diamond 
atoms more accurately than the ZBL potential used in 
CTRIM. This leads to a higher number of channeled 
projectiles. 



B. Annealing 
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Figure 3: Depth distribution of single implanted nitrogen 
atoms with an incident angle of 0° and with an implantation 
energy of 4keV; red: MD simulation, grey: CTRIM simula- 
tion. The black and blue curves are Gaussian fits summarized 
in a cumulative curve (green). 



A broad distribution can be extracted from simulation 
where a "not channeled" and a "channeled" part of 
the implanted nitrogens appears. In the MD case, 
nitrogens are implanted at a depth of — 460A. The 
data can be fitted by using two Gaussian peaks at 
d c i = 75 A (first peak: black curve) and d C 2 = 2 10 A 
(second peak: blue curve) . The first Gaussian peak 
has a FWHM of 77 A whereas the second Gaussian 
peak has a FWHM of 177 A. The number of implanted 
nitrogen atoms at d c ± is 21 and at d C 2 10. In addition 
the simulations show that 37% of implanted nitrogens 
ending up on a lattice site and 63% end up as interstitials. 

The CTRIM results (figure [3|) differ from the MD 
outcome. Here the number of projectiles was chosen 
to be 10 6 ; the CTRIM results are scaled in order to 
make them comparable with the MD results. The 
"not channeled" part shows a 25% narrower and the 
"channeled" part shows a 33% broader distribution 
than in the MD simulation. The maximum of the "not 
channeled" peak is shifted by 25A (33%) towards lower 
values whereas the maximum of the "channeled" peak is 
shifted by 20 A (10%) towards higher values compared to 
MD results. This might be due to the fact that during 
the implantation a higher number of nitrogens enter 



After the implantation step it is unlikely to find a next 
neighbor pair of a nitrogen atom and a vacancy. Not 
until the annealing step, when vacancies become mobile 
and diffuse through the lattice, can such pairs be formed. 
Several studies on diffusion near to the surface have been 
done in recent years22r— . Here we modeled the diffusion 
process during annealing and investigated its influence 
on the NV center yield. 

The statistical probability for vacancies, created during 
the implantation, to migrate to a next neighbor position 
of the implanted nitrogen was examined by Monte Carlo 
simulations. Our Analysis included the most important 
loss mechanisms for vacancies, such as migration to the 
surface or recombination with interstitials. The processes 
in which a self-interstitial atom replaces an on-site nitro- 
gen at a lattice site and turns the latter into an interstitial 
nitrogen atom 3 - 3 - was also included. In contrast to on-site 
nitrogen atoms which hardly move at HOOK, interstitial 
nitrogen atoms are mobile during annealin g 33 ! 34 . The 
initial distribution after implantation was assumed to be 
40% on-site nitrogen and 60% between lattice sites, as 
suggested by the MD simulations discussed above. Fig. 
H] shows the percentage of nitrogen atoms that are found 
on lattice sites after the annealing. 

The number increased slightly and is the same for all 
energies within error bars. Hence only about half of the 
implanted nitrogen atoms are available for NV center cre- 
ation in the first place. This is in good accordance with 
the maximum yield found by Pezzagna et al£ for MeV 
implantations. It suggests that for high energies every 
implanted nitrogen that finally ends up on a lattice site 
is turned in to a NV - center. Figure [5] shows the per- 
centage of initially created vacancies that are lost at the 
surface or due to recombination. By considering intial 
SRIM input (5-50keV), as expected, at shallow implan- 
tations more vacancies (48%) are lost at the surface than 
for deep ones (26%). This is because of the short distance 
to the surface that makes it more likely for vacancies and 
self-interstitials to randomly migrate to the surface where 
they are lost, instead of recombining. Except for very low 
implantation energies (<20keV) the loss due to recombi- 
nation with interstitial carbon is significantly higher. It 
increases from 44% at 5keV to 62% at 50keV while the 
loss at the surface decreases from 49% to 26%. The line 
crossing of the vacancies lost at the surface and due to 
recombination define a threshold. For implantation ener- 
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Figure 4: Number of on-site nitrogens after annealing depend- 
ing on the implantation energy. The inset shows the time 
evolution of on-site nitrogens during annealing after 50 keV 
implantations with an initial 40% on-site and 60% interstitial 
distribution. The number of on-site nitrogen atoms increases 
already within the first simulation steps which are not shown 
in the inset. 
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than the results based on SRIM data. IMD seems to yield 
defect structures with smaller distances between vacan- 
cies and the corresponding interstitials which therefore 
are more likely to recombine. The inset in figure [5] shows 
the percentage of remaining vacancies that are still avail- 
able for NV center creation. It can be seen that at the 
end of the simulation the vacancy diffusion within the 
crystal has almost completely come to an end. On aver- 
age 90% of the created vacancies during the implantation 
process are lost after 2 hours of annealing. 
Figure [6] shows the yield of NV centers of any charge 
state that is to be expected statistically. It is about the 
same for all investigated implantation energies, roughly 
25% of implanted ions. 




Figure 6: Number of NV centers created after implantation 
and annealing. 



IV. SUMMARY AND DISCUSSION 



Figure 5: Number of vacancies which are lost at the surface 
(red) or lost due to recombination processes (black) after 2 
hours of simulated annealing. The inset shows the number of 
vacancies that are still available for NV center creation over 
time. 



gies lower than 7keV the energy loss is dominated by the 
surface process. For higher implantation energies recom- 
bination process reduce the available vacancies to form 
NV centers. The simulations show that the recombina- 
tion of vacancies and interstitials happens within the first 
simulated seconds of annealing because of the high mo- 
bility of carbon interstitials. Hence the number of vacan- 
cies that are left to possibly migrate towards a nitrogen 
atom is less than half the number that is suggested by 
implantation simulations such as SRIM, which does not 
take recombination into account. The results at 4keV 
implantation energy based on IMD implantation simula- 
tions show an even higher probability of recombination 



In this paper we present results using MD, DFT and 
MC simulation techniques to investigate the NV center 
yield in diamond. We could show that during nitrogen 
implantations into diamond even at low implantation 
energies (4keV) channeling effects occur and have to be 
taken into account for accurate placement in the carbon 
matrix. Furthermore, the annealing of the implanted 
diamond structures lead to a final NV center yield of 
about 25%. This result is significantly higher than the 
NV - yield seen experimentally (0.01 — 0.1%)^ with 
implantation energies <10keV. The discrepancy could 
be explained in the final charge state of the shallow 
implanted NV centers. Experimentally only the NV - 
and the NV° can be detected. Surface effects or damage 
structures in the vicinity of the NV center could change 
the charge of the NV - ; leaving the optically inactive 
NV + state. Our simulations do not consider charges and 
therefore give an average value for NV center creation 
of all charges. Introducing additional electrons to the 
system which can be captured by the NV centers to 
form the NV - may potentially increase the very low 
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experimental yield. 

Further investigations on the implantation and the 
annealing process are necessary. The IMD simulation 
package considers only the nuclear stopping a valid 
approximation for the low energy implantations used, 
whereas the electronic stopping is not taken into account. 
In addition the NV yield depending on the annealing 
temperature, the intrinsic nitrogen content and the im- 
plantation densities have yet to be studied. In particular 
efficient simulations for the implantation with higher 
energies require simulations of larger systems which is 



out of reach with current codes and computing powers. 
Adding precise calucations of the electronic states of 
defect structures would yield a comprehensive picture of 
the defect center and its surrounding necessary for the 
understanding of the generation process of larger spin 
arrays. 
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